clc;clear;clearvars;
% 随机生成20个数据
num_initial = 20;
num_vari = 60;
% 搜索区间
upper_bound = 5.12;
lower_bound = -5.12;
% 迭代时计算y值总的个数
eval_num = 12000;
% K表示划分子空间的个数,sub_num表示每个子空间的维度
K = 5;
sub_num = num_vari / K;
% 算法的迭代次数
iter = eval_num / (2 * K);
w = 1;
% 随机生成20个数据,并获得其评估值
sample_x = lhsdesign(num_initial, num_vari).*(upper_bound - lower_bound) + lower_bound.*ones(num_initial, num_vari);
sample_y = Rastrigin(sample_x);
Fmin = zeros(iter, 1);
aver_Fmin = zeros(iter, 1);

for n = 1 : 100
    n1 = 1;
    % 初始化一些参数
    pbestx = sample_x;
    pbesty = sample_y;
    % 当前位置信息presentx
    presentx = lhsdesign(num_initial, num_vari).*(upper_bound - lower_bound) + lower_bound.*ones(num_initial, num_vari);
    vx = sample_x;
    [fmin, gbest] = min(pbesty);
    global_best_x = pbestx(gbest, :);
    fprintf("n: %.4f\n", n);
    % fprintf("iter 0 fmin: %.4f\n", fmin);
    for i = 1 : iter
        index = randperm(num_vari);
        for i1 = 1 : K
            % 方法1：随机分组
            ind = index((1 + (i1 - 1) * sub_num) : i1 * sub_num);
            % ind = ((1 + (i1 - 1) * sub_num) : i1 * sub_num);
            r = rand(num_initial, sub_num);
            % pso更新下一步的位置,这里可以设置一下超过搜索范围的就设置为边界
            vx(:, ind) = w.*vx(:, ind) + 2 * r .* (pbestx(:, ind) - presentx(:, ind)) + 2 * r .* (pbestx(gbest, ind) - presentx(:, ind));
            vx1 = vx(:, ind);
            vx1(vx1 > upper_bound) = upper_bound;
            vx1(vx1 < lower_bound) = lower_bound;
            vx(:, ind) = vx1;
            presentx(:, ind) = presentx(:, ind) + vx1;
            presentx1 = presentx(:, ind);
            presentx1(presentx1 > upper_bound) = upper_bound;
            presentx1(presentx1 < lower_bound) = lower_bound;
            presentx(:, ind) = presentx1;

            presentx2 = repmat(global_best_x, num_initial, 1);
            presentx2(:, ind) = presentx1;
            presenty = Rastrigin(presentx2);

            % 更新每个单独个体最佳位置
            pbestx1 = pbestx(:, ind);
            pbestx1(presenty < pbesty, :) = presentx1(presenty < pbesty, :);
            pbestx(:, ind) = pbestx1;
            pbesty(presenty < pbesty, :) = presenty(presenty < pbesty, :);
            
            % 更新所有个体最佳位置
            [fmin, gbest] = min(pbesty);
            global_best_x = pbestx(gbest, :);
            %disp(fmin);
        end
        % 方法2：自适应权重,种群个数为5
        weight_pbestx = ones(5, sub_num);
        weight_vx = ones(5, sub_num);
        weight_pbesty = repmat(fmin, 5, 1);
        [weight_fmin, weight_gbest] = min(weight_pbesty);
        weight_presentx = rand(5, sub_num);
        for i2 = 1 : 20
            r3 = rand(5, sub_num);
            global_best_x1 = repmat(global_best_x, 5, 1);
            weight_vx = 0.72 * w.*weight_vx + 2 * r3 .* (weight_pbestx - weight_presentx) + 2 * r3 .* (weight_pbestx(weight_gbest, :) - weight_presentx);
            weight_presentx = weight_presentx + weight_vx;
            for i3 = 1 : K
                global_best_x1(:, (1 + sub_num * (i3 - 1)) : sub_num * i3) = weight_presentx(:, i3) .* global_best_x1(:, (1 + sub_num * (i3 - 1)) : sub_num * i3);
            end
            global_best_x1(global_best_x1 > upper_bound) = upper_bound;
            global_best_x1(global_best_x1 < lower_bound) = lower_bound;
            weight_presenty = Rastrigin(global_best_x1);
            weight_pbestx(weight_presenty < weight_pbesty, :) = weight_presentx(weight_presenty < weight_pbesty, :);
            weight_pbesty(weight_presenty < weight_pbesty, :) = weight_presenty(weight_presenty < weight_pbesty, :);
            % 更新所有个体最佳位置
            [weight_fmin, weight_gbest] = min(weight_pbesty);
            %disp(weight_fmin);
        end
        fmin = weight_fmin;
        for i4 = 1 : K
            global_best_x(1, (1 + sub_num * (i4 - 1)) : sub_num * i4) = weight_pbestx(weight_gbest, i4) .* global_best_x(1, (1 + sub_num * (i4 - 1)) : sub_num * i4);
        end
        global_best_x(global_best_x > upper_bound) = upper_bound;
        global_best_x(global_best_x < lower_bound) = lower_bound;
        pbestx(gbest, :) = global_best_x;
        pbesty(gbest, :) = fmin;
        Fmin(n1, 1) = fmin;
        n1 = n1 +1;
        %fprintf("iter %d fmin: %.4f\n", i, fmin);
    end
    aver_Fmin = aver_Fmin + Fmin;
end
aver_Fmin = aver_Fmin ./ 100;
plot(aver_Fmin);